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In the companion paper of this set (Capitan and Cuesta 20101 we have developed a full analytical treatment of the model of 
species assembly introduced in 'Capitan et al. ( 2009[ l. This model is based on the construction of an assembly graph containing 
all viable configurations of the community, and the definition of a Markov chain whose transitions are the transformations of 
communities by new species invasions. In the present paper we provide an exhaustive numerical analysis of the model, describing 
the average time to the recurrent state, the statistics of avalanches, and the dependence of the results on the amount of available 
resource. Our results are based on the fact that the Markov chain provides an asymptotic probability distribution for the recurrent 
states, which can be used to obtain averages of observables as well as the time variation of these magnitudes during succession, 
in an exact manner Since the absorption times into the recurrent set are found to be comparable to the size of the system, the 
end state is quickly reached (in units of the invasion time). Thus, the final ecosystem can be regarded as a fluctuating complex 
system where species are continually replaced by newcomers without ever leaving the set of recurrent patterns. The assembly graph 
is dominated by pathways in which most invasions are accepted, triggering small extinction avalanches. Through the assembly 
process, communities become less resilient (e.g., have a higher return time to equilibrium) but become more robust in terms of 
resistance against new invasions. 
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^1. Introduction 

Understanding the mechanisms leading to species assembly 
in ecological communities is a challenging issue. In particular, 
assembly models have been used to understand the observation 



that natural communities are both complex and stable ( McCann 
[20001 |Dunne||2006| . 

Assembly models try to mimic the sequential arriving of 
rare species (invaders) to which natural communities are sub- 
jected. Standard assembly models ( |Drake|[T990[|Law and Mor^ 



ton 



1993[|1996) use "species pools" as (finite) sets of potential 



invaders. Pools are usually defined by labeling species accord- 
ing to some niche variable (usually a species trait like body 
size) and then drawing randomly their interactions from pre- 
determined probability distributions ( [Law and Mortonj |1996| l. 
Sequential invaders of any given resident community are se- 
lected from the pool at each invasion attempt, and the resulting 
community after the invasion can be determined according to 
some population dynamics. For models using Lotka-Volterra 



equations, the permanence (Hofbauer and Sigmund 1998 1 of 
the invaded community is a suitable criterion which determines 
the same final community as the numerical integration of the 
equations (Morton et al. 1996| l. 
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The most remarkable results of previous assembly models 
are: (i) a final end state is eventually reached, which can be 
either a single community or a cycle involving several com- 



munities ( {Morton and Law 1997 1, (ii) average species richness 



(complexity) increases with successional time ( |Post and Pimm 



[1983; Drake 1990; |Law and Morton| [T996l l, and (iii) stability, 
understood as resistance against invasions, also increases with 
time ( |Case| |1990[ |Law and Morton] |1996[ [Morton and Law[ 



1997) 1. Thus assembly models conform a well-founded theo- 



retical framework that provides a positive relationship between 
stability and complexity in model communities. 



In a previous paper (Capit an et al.||20()9] l we have provided a 
picture of the assembly process of an ecosystem as a Markov 
chain evolving in a certain configuration space. This space 
is made of all viable communities for a given set of parame- 
ters (resource saturation, interspecific competition, consump- 
tion rates. . . ). The invasion process by a new species induces 
transitions as a result of the perturbations created in the com- 
munity by the newcomer. The process drives the community to 
an end state resistant to invasions. For some parameter values 
this end state is just a single uninvadable community. For the 
remaining values, the end state is formed by a set of communi- 
ties with equal number of trophic levels and similar number of 
species per level, which transform into each other as a result of 
new invasions. In this set, communities can always be invaded 
but they never abandon the set. These complex end states are 
a generalization of the end cycles found in previous assembly 
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models (Morton and Law 1997| l, and the fact that they had not 
been observed so far is probably due to limitations in the pool 
of invaders of these previous models. In the preceding paper 



( Capitan and Cuesta 2010j we have shown that the existence of 
these complex end states is a result of a top predator attempting 
to invade a community when its establishment is not allowed by 
the parameters of the model. 

Our model recovers the main findings of previous assembly 
models (Post and Pimm! |T983) |Drake| [T9901 [Case] [T9901 |Law| 
[and Morton, .1993, .1996l l, such as the resistance of end states 
against invasions, or the increase of complexity (biodiversity) 
along the assembly (Capitan et al. 2009 | l. Its main virtue is 



therefore being sufficiently simple so as to allow mapping out 
all assembly pathways, thus providing a global picture of the 
assembly process. 

Despite recovering the above similar results, there are impor- 
tant differences between our work and early models. First, the 



niche variable in our model is simply the trophic level ( Capitan 



et al. 2009 j , which renders our species pool infinite (in con- 



trast to most previous models; but see Post and Pimm ( 1983 1 for 
an exception). However, interactions in the pool are averaged 
over each trophic level under a species symmetry assumption 
(Capitan and Cuesta 2010[ l, which decreases substantially the 
number of different assembly pathways. Second, in our model 
the permanence of the final community is guaranteed because 
we are able to show that equilibrium communities are globally 
stable ( [Hofbauer and Sigmundj [T998| l under the assumption of 
neutraUty within each trophic level. And third, in standard as- 
sembly models magnitudes are averaged over a set of stochastic 
realizations of the process of sequential invasions, where in- 
vaders are randomly chosen from the species in the pool not yet 
present in the community. Since we are able to map out all the 
invasion pathways for this model, we do not need to resort to 
average magnitudes over realizations but we can calculate them 
exactly (Capitan et al. 2009 ) l. Even for our simple model, the 
number of possible pathways is too high to be accounted by 



through simulation [see details in Capitan et al. ( 2009) ]. This is 
one of the main advantages of our model with respect to former 
ones, which in turn allows us to establish its independence on 
history. The uniqueness of the end state for these kind of mod- 
els was not proven until now, although it was already found that 
most of the simulated assembly sequences led to a single set of 
final communities ( Morton and Law 1997| l. 



In the first paper of this suite (Capitan and Cuesta 2010 1 we 



have performed a detailed analysis of the analytical properties 
of the Lotka-Volterra population dynamics underlying our as- 
sembly model, as well as the stability properties of the interior 
equilibria from the dynamic point of view. Our communities 
represent a mean-field version of trophic networks: the feeding 
relations are assumed to take place only between contiguous 
trophic levels and the strength of each interaction is averaged 
to a uniform value. This assumption of symmetry allowed us to 
simplify the differential equations, showing that in our model 
the set of species numbers at each level is enough to 

determine the equilibrium densities and the dynamics of a com- 
munity with L trophic levels. 

The present paper presents new results of the assembly pro- 



cess ranging from structural properties of the building of stable 
communities to dynamical properties of the stochastic invasion 
process. The transition matrix associated to the chain allows 
us to define a stationary distribution of probabilities over the 
assembly graph. We will use that distribution to calculate aver- 
ages of biologically relevant observables in the ecosystem, like 
the average number of species, the total population density, etc., 
and even to obtain distributions of certain magnitudes like the 
fraction of extinct species after a destructive invasion. But this 
is not the only advantage of our approach. The transition ma- 
trix also provides the time evolution of the probability, so the 
dependence in time of any magnitude can be obtained exactly. 
Due to its simplifications, the model reduces the number of pos- 
sible communities to a finite graph. Once we have that graph, 
the theory of finite Markov chains can be exploited to obtain the 
dynamical properties of the assembly process. 

The paper has been organized as follows. In Section [2] we 
revisit the definition of the Markov chain associated to the as- 
sembly process and show that all communities can be classified 
into transient or recurrent. Section [3] is devoted to discuss the 
main statistical results that can be obtained for our model, for 
instance, the asymptotic distribution within the complex end 
states (Section 3.1 1, the dependence of averages upon varia- 



tion of the resource saturation (Section 3.2 1, the dependence of 



the results upon variation of the parameters of the model (Sec- 
tion |3J]l, the average time to reach the end state (Section [3 .4| l, 
the statistics of avalanches of extinctions caused by invasions 



(Section 3.5 i or the variation of biologically relevant averages 
with successional time (Section |3.6|l. We finally discuss our 



findings and their implications in Section|4] 

2. The assembly process as a Markov chain 

Let us start with a detailed description of the Markov chain 
associated to our assembly model In the first paper of this suite 
(Capitan and Cuesta 2010| l we showed that, solving the linear 
system that defines the interior equilibrium point of the Lotka- 
Volterra equations (i.e., the system obtained by equating to zero 
the per-capita population growth rates «, /«,, with / running over 
the set of all species) we can determine all viable communities 
(i.e., those whose population densities are above a certain ex- 
tinction threshold tic > 0) that are compatible with a given set 
of parameters. Although in principle the population model al- 
lows for infinitely many species at each level, it turns out that 
the set of viable communities is finite. This is a consequence of 
the existence of the extinction threshold, that precludes the ap- 
pearance of arbitrarily small populations (an unrealistic feature 
of deterministic population dynamics models). There is another 
limitation due to the finite amount of abiotic resource that main- 
tains our model communities. In our previous work we modeled 
this resource with a linear functional response with a saturation 
value R (Capitan et al.| |2009[ [Capitan and Cuesta| [2010 ). R 
accounts for the amount of resource that would be reached in 
the absence of consumers. On the one hand there is a maxi- 
mum number of levels allowed for a given resource saturation 
R (Capitan and Cuesta 2010[ l; on the other hand population 
densities at each level decrease as 1 / S(, with S( the number of 
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Figure 1 : (Color online) Assembly graph obtained for a value of the resource 
saturation R = 80. It is made of 39 communities (nodes), each of them with 
either one or two trophic levels. Transitions shown with a black arrow indicate 
that the invasion is accepted, and those with a red arrow refer to a reaiTangement 
in the resulting community after the invasion. Transient nodes are filled in red, 
and recurrent nodes are filled in dark blue. In this case, the final end state 
of recurrent communities comprises 3 communities forming an end cycle like 
those found by |Moi1on and Law]jl997) . Labels of each node show the species 
numbers {.51,^2} of each trophic level in the community. (For the choice of 
parameters see Section[X2]) 



species in that level ( Capitan and Cuesta 2010| l, so we can have 
populations infinitely close to zero. Therefore, the existence of 
the extinction threshold renders the set of communities under 
consideration finite, and then the associated Markov chain has 
a finite number of states. Besides this being a more realistic 
description of an ecosystem, it also drastically simplifies the 
analysis of the assembly process. 

Thus for any choices of parameters there is a finite set of 
viable communities — that we denote by ;F. There will be a 
link from community / to community j of the set provided the 
former is transformed into the latter as a result of an invasion. 
Invasions are assumed to occur at a uniform rate f . We assume 
that the typical dynamical time is much smaller than the 
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Figure 2: (Color online) Same as Figure [T]for R = 50. The total number of 
communities in this graph is 16, and they have up to 2 trophic levels. An end 
state with a single community (with occupancies {4, 2}) is reached in this case. 



mean time between invasions, so that communities are always 
at equilibrium when an invasion occurs [for the validity of this 
assumption see the discussion in Capitan and Cuesta ( 2010| l] 



The population of the invader is assumed as small as possible, 
i.e. equal to ( Roughgarden] [T9741 |Tureni| [TMT] ). 

Consider a community i G 'F', with L trophic levels, at 
its equilibrium point. Potential invaders are species of level 
{ - 1 , . . . , L + 1 (species of higher levels would not be able to 
feed from the existing levels). We randomly choose { and intro- 
duce a new species at level £ of the community /. The extended 
community will evolve to the interior equilibrium correspond- 



ing to the new number of species if + 1 at level { ( Capitan and 



Cuesta 2010 1. If this equilibrium is viable, then the invader is 



accepted in the resident community. The new community j will 
also be in 'F and a directed link will go from / to j. If the new 
equilibrium is not viable then we apply the procedure discussed 
in jCapitan and Cuesta ( 2010 ) to determine what is the sequence 
of extinctions until the community becomes viable. Two things 
can happen: either the first level to fall below the extinction 
threshold is the invaded one, or it is another one. In the former 
case the invader is simply rejected and the community remains 
unaltered; in the latter, the extinction sequence will lead to a 
new community k eT, and a link will go from / to k. 

The assembly graph, Q, is defined as the connected compo- 
nent containing the empty community, 0, of the directed graph 
whose nodes are the elements of T and whose links are the 
transitions obtained by the invasion process just described. Ob- 
viously, the way to construct Q is to start off from 0, and pro- 
ceed by attempting all possible invasions for every community 
reached along the assembly process [see Figures [T] and [2] for a 
graphic representation of simple assembly graphs, for a larger 



value of R see Figure 1 in Capitan et al. (2009) ]. The exhaustive 
characterization of the set of nodes in ^ is a bit demanding. For 
example, the storage of all the communities appearing along 
the process has been carried out by using a binary tree (Knuthj 
|1997| l, by exploiting the mapping between any configuration 
{lyj^j and a binary number (we simply concatenate the binary 
representations of each species number S[ to form a branch of 
the binary tree). Despite this, we have been able to analyze 
graphs with around 10^ communities within. 
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Figure 3: (Color online) Graphs of the recurrent subset for values of the resource saturation R = 990 (a) and R = 390 (b) (other parameters are set as in Figures[T]and 
|2j. The first graph contains 10 communities with 4 trophic levels, whereas the second has 30 communities with 3 levels. For a graphical representation of a recun'ent 
subset more complex than these see |Capitan et al.H2009j l. All the communities in both cases have three trophic levels. The diameter of the nodes is proportional to 
its asymptotic probability. Black arrows show accepted invasions and yellow ones those causing a reconfiguration (the thickness of each line is proportional to the 
relative number of extinct species). Labels indicate the number of species in each trophic level. 



The connection of the species assembly process with a 
Markov chain on the graph § amounts to assigning certain tran- 
sition probabihties to each Hnk of the assembly graph. We de- 
fine these probabilities in a simple way. Invaders arrive at each 
community at a constant rate ^, independent of the level of in- 
vasion, and the stochastic process is updated in discrete time 
(each time unit is the average time elapsed between consecu- 
tive invasions). Thus, if / and j are two nodes of connected 
by a link, we assign it the transition probability 



(1) 



where 6ij = 1 if / = j and otherwise. The matrix elements Qjj 
are given by 



(2) 



where is the number of different invasions of i that lead to j 
and L -H 1 is the number of different invasions of /, provided it 
has L trophic levels. Therefore, the probability of the transition 
i j between different communities is proportional to the rel- 
ative frequency of the transition among all the possible transi- 
tions starting from /, the invasion rate being the proportionality 
constant. The diagonal of Q is chosen so that P - {Pij) is a 
stochastic matrix. 

Since the diagonal elements of the transition matrix P are 
non-zero, the Markov chain can not be periodic ( |Karhn and] 
Taylor 1975 1. As the set of viable ecosystems T is finite, P de- 
fines the transition matrix of a finite, aperiodic, Markov chain. 
The states of one such chain are either transient or recurrent 
(Karlin and Taylor 1975[ ). There can be one or several subsets 



of recurrent states, the chain being ergodic in each of them. Ev- 
ery recurrent subset is a different end state of the assembly pro- 
cess. The end state of an ecosystem will be history-dependent 
only if there are at least two such recurrent subsets. Ergodic- 
ity implies that there is a stationary probability distribution on 
the states of these subsets which determines the frequency with 
which the process visits each of them. [For a full account on 



Markov chains see e.g. Karlin and Taylor ( 1975])]. Our model 
only exhibits a unique recurrent set for any given set of param- 
eter values ( [Capitan et al. 2009 1. 

This concludes the definition of the Markov model for the 
assembly process. As we are able to compute the whole tran- 
sition matrix P, we have a complete and exact characterization 
of the assembly process. In particular, by selecting an initial 
state for the Markov chain (in our case the process always starts 
off from 0), we can obtain the evolution of any magnitude — 
numerically but exactly — without resorting to taking averages 
over realizations of the process. In the following section we 
will discuss in detail the results that can be obtained from the 
analysis of the Markov chain [a brief account of which were 
reported in Capitan et al. ( 2009| l]. 



3. Results 

3.1. Asymptotic distribution 

To separate transient and recurrent states, we have applied 
an algorithm provided by 'Xie and Beerel ( 1998'). Notice that 
the characterization of transient and recurrent states in a finite 
chain depends only on the graph, not on the transition proba- 
bilities. Only one subset of recurrent states was found for each 
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parameter value interpretation 

R 10 < R < 1700 Saturation value, in the absence of predation, of the abiotic resource abundance 

a 1 Average mortality rate of consumers 

y_ 5 Average rate of decrease in preys population caused by their being predated 

y+ 0.5 Average rate of increase in predators population due to feeding 

p < p < 1 Relative magnitude between intra- and interspecific competition 

He 1 Extinction threshold 



Table 1 : Summary of parameters of the model and ecological meaning of each one of them. 
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Figure 4: Distribution of asymptotic probabilities n within each recurrent set, 
for several values of R. Communities are labelled in decreasing order of prob- 
ability. The distributions look exponential with a cutoff. 



set of parameters. Let denote the subgraph of formed by 
this ergodic set. Figure [3] shows two examples of these sub- 
graphs. The particular transition probabilities assigned to each 
link would determine the asymptotic probability distribution 
within the recurrent set, but not the subset of nodes contained 
in it. 

In order to calculate the asymptotic probability nt for a com- 



munity i G we need to solve the linear system n - nP (Karlin 



and Taylor 1975 1, in other words, the vector n of probabilities 
is the left eigenvector of the matrix P with eigenvalue 1 . Since 
our graphs are very sparse, standard numerical techniques for 
solving sparse systems have been applied. The eigenvector is 
normalized to satisfy the condition YiieK^i - 1- Obviously, we 
only need to solve this system for the subgraph corresponding 
to the recurrent set, since by definition the asymptotic proba- 
bility TTi - for any transient state /. Note that our matrix of 
transition probabilities ([TJ reduces the condition to be satisfied 
by n to jiQ - 0, i.e., tt is a left eigenvector of Q with eigenvalue 
0. It is worth noticing that neither the asymptotic distribution, 
nor the recurrent subset depends on the invasion rate. 

We can thus obtain a probability distribution for each recur- 
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Figure 5: (Color online) Total number of communities (black circles, below) 
and transitions (red crosses, above) in the Markov chain as a function of the 
resource saturation R. 



rent set. Having this probability distribution is therefore equiv- 
alent to defining a statistical mechanics over the set of viable 
communities, if we regard ;F as the phase space of our system. 
In Figure [4] we have plotted the histogram of probabilities for 
several values of the resource saturation R (for these values the 
number of communities in each set is larger than 10^). Com- 
munities are labelled in decreasing order of probability. These 
distributions are found to be roughly exponential over several 
orders of magnitude, this meaning that only a small number 
of communities (in general very similar to each other) occur 
with high probability. These are the communities in which it 
is more likely to find the ecosystem. Nonetheless ergodicity 
implies that all communities in the end state are visited with 
non-zero probability. The ecosystem is thus in a complex state, 
with fluctuating species numbers in each level due to some in- 
vasions being accepted and some others causing avalanches of 
extinctions. 

This distribution can be used to calculate the asymptotic aver- 
age over of any relevant magnitude M, defined for every com- 
munity, like for instance the average number of species, the total 
population, etc. We just need to evaluate {M}ii = TjieR ^iMj. 

3.2. Dependence with the resource saturation 

All results presented here have been obtained with a death 
rate a - 1, an extinction threshold = 1, an average rate of 
increase in predators population per predation event 7+ = 0.5, 
and an average decrease in reproduction rate of preys per preda- 
tion event j- - 5 (see Capitan and Cuesta ( 2010] l for details on 
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Figure 6: Number of viable states that are not reachable trough invasions start- 
ing from the empty community 0, AN = Nf- - Ng. Typically there is an accu- 
mulation of these communities in the regions where R allows a top predator in 
them. 



Figure 7: (Color online) Statistics of links con'esponding to accepted invasions 
(red squares), rejected invasions (blue crosses) and invasions that lead to a re- 
configured community through an avalanche of extinctions (black circles). Dot- 
ted lines correspond to the onsets of acceptance of a new trophic level, and 
dashed lines to the beginning of the regions of complex end states. 



the way these parameters enter the population equations of the 
model and Table [T] for a brief summary of them). The assump- 
tion of 7+ <K y_ is ecologically sound, because many preys 
must be consumed to produce a new predator, while loosing one 
prey requires a single predation event. A common choice for 



the energy transfer between trophic levels is about 10% (Pimm 



|1991) l. We have checked that the model is robust against varia- 
tions of these parameters within reasonable bounds. 

In most cases we have taken the ratio of direct inter- to in- 
traspecific competition p = 0.3. However, we have explored 



the effect of this parameter in detail in Section 3.3 



We have obtained all assembly graphs in a range of resource 
saturations that goes from R = 10 up to 7? = 1700 with incre- 
ments AR = 5. No viable community is found below R - 10. 
The number of communities Ng in these graphs goes from just 
one (for R = 10) up to about 10^. We have found empirically 
that both this number and the total number of transitions in each 
graph grow roughly as e"^, see Figure [s] The maximum num- 
ber of trophic levels that are allowed up to /? = 1700 is 5. 

We have checked whether the set of communities in the as- 
sembly graph is the whole set Given the estimation of the 
resource values that allow a maximum number of levels Lmax 



[see Capitan and Cuesta (2010|], we have checked the viability 



of all possible combinations of species numbers {s(}^2'i with 
imax + 1 levels up to a total number of species 5 max equal to 
twice the maximum number of species allowed for that R value. 
Since there is a huge number of these combinations when R in- 
creases, we have checked this up to /? = 700. Figure |6] shows 
the difference AA^ - Nf- - Ng. In nearly all cases the set of 
communities in the assembly graph is but we have found 
several instances — all of them near the values of R at which 
new levels arise — in which contains communities not reach- 
able through the assembly process, just like in the experiment of 
Warren et al. (2003). The largest difference is found at 7? = 470, 
where Ng - 4800 and AA^ = 375, so the highest relative differ- 
ence reaches 8%. 

For each R we determined the number of recurrent states of 
the chain (see Figure 3 in Capitan et al. ( 2009 [ l for a plot of this 



number as a function of R). We always found a single connected 
graph, which implies that the end state of the assembly process 



does not depend on history for this model (^Drake 1990 1. This 
is consistent with the results of Morton and Law "( |1997 1 as well 
as the experiments of [Warren et al. ( 2003| l. There are values 
of R for which this set consists of a unique absorbing state (or 
just a few, sometimes forming a cycle), but when R is reaching 
the values at which a new trophic level appears, the size of this 
set increases considerably (the largest set found contains around 
1800 communities; a tiny fraction of the whole assembly graph, 
anyway). After crossing these values the size of the recurrent 



set drops again down to just one absorbing state. Morton and 



Law ( 1997 1 also obtained complex end states in 6 out of the 



80 pools they explored, with a number of communities ranging 
from 6 to 138. 

In Figure [7] we show the fractions of links in the assembly 
graph corresponding to invasions that are accepted, rejected, 
or cause a reconfiguration in the system through a sequence of 
extinctions. The most frequent case is the acceptance of the 
invader, although there are around 20% of rejections and recon- 
figurations. We can see an increasing trend to reconfigurations 
when R corresponds to a complex end state [see Figure 3 in 
Capitan et al. ( 2009[l]. T he invasibility criterion discussed in 



Capitan and Cuesta ( 2010 1 explains why we observe an increas- 



ing number of rearranged communities in these regions. 

As for the dynamic stability (resilience), we can measure the 
return time, i.e., the mean time that a perturbed ecosystem needs 
to restore equilibrium ( Pimm and Lawton 1977| l, averaged over 
the probability distribution of the stationary state. It can be ob- 
tained as Tr = "^m'ax' whcrc /Imax IS the largest real part of 
the eigenvalues of the linear stability matrix — which is always 
negative in our communities since they are globally stable. We 
observe that this time is roughly independent on the end state, 
being approximately constant as a function of the resource sat- 
uration R (see Figure[8^). 

For each end state, regardless on whether it is an absorbing 
community or a recurrent set, we have calculated another aver- 
age. In FigurelSb we show the dependence of the total popula- 
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Figure 8: (a) Mean return time in tlie stationary state vs. resource saturation R. 
The behavior is approximately constant, except for the region of low resources, 
where the graphs contain less communities and there is more variability, (b) 
Mean population density of the community vs. R. 



Figure 9: Number of communities in the recurrent sets vs. the resource satura- 
tion varying direct competition (upper panel, p = 0; middle panel = 0.3; lower 
panel, p = 0.7). Observe the decrease in Nn as competition increases, and the 
increase of the values of R at which a new level arises. 



tion of a community, B - Yti^i stp^, averaged over the recur- 
rent set as a function of the resource saturation R (p^ denotes 
the equihbrium density of the species at level f). The depen- 
dence is practically linear, except for some dips near the onset 
of emergence of a new level. 

3.3. Dependence on the parameters 

We have already mentioned that the model results are not 
qualitatively influenced by variations of its parameters. For ex- 
ample, we have studied the model dependence with respect to 
direct competition (Figure |9]l. In the absence of interspecific 
competition (p = 0) levels are filled up more easily, so the num- 
ber of communities in the recurrent set increases with respect to 
the results reported so far The effect of increasing direct com- 
petition is to reduce the number of ecosystems in these sets, and 
to increase the resistance to the appearance of a new level in the 
end state for the same values of resource saturation. Thus, the 
global behavior of the number of communities as a function of 
R turns out to be similar, up to scale factors, to that obtained in 



Capitan et ai:] ( |2()()9j ). Figure 3. 

The particular case p - \ (interspecific competition equal to 
intraspecific competition) is qualitatively different. Fixing p - 
1 transforms the community into a trophic chain. All species 
can be grouped into a single one with population A^'^ = s^n^ 



[i.e. th e Lotk a-Volterra equations of this system ( Capitan and 
Cuesta 2010 1 are closed in the variables A^^]. But the impli- 



cations of this assumption are stronger. Even if the distinction 
between species becomes meaningless, one can formally keep 
the identities and treat them as different. But then it is easy to 



show that any invasion attempted at a level akeady occupied by 
at least one species will be unsuccessful because the population 
of the invader ends up below n^. In fact, according to Eq. (12) 



of Capitan and Cuesta (2010 1, the initial per capita growth rate 
of an invader at level £ is —n^ and the equations for level I and 
for the invader coincide. Hence jN^ - h/n, n being the abun- 
dance of the invader. This asymptotically yields 



P = 



A^^(O) 



(3) 



since n(0) = «£. (p and are the invader and the total ^-level 
densities at equilibrium after the invasion, respectively). Now, 



the linear system (13) in Capitan and Cuesta (20101 for the in- 



terior equilibrium point before the invasion is exactly the same 
after the invasion replacing A^^(O) hy p + P^. This fact, together 
with ([3|, yields p - ncN^(0)/[nc + A^^(O)] < n^. Since the pop- 
ulation of the invader initially decreases, according to our ex- 



tinction procedure (Capitan and Cuesta 2010 1 the invader goes 
extinct. 

Thus the assembly graph § becomes trivial. Using the nota- 
tion {sfl^^i for each community, Q is simply 



|1,0,...,0) ^ {1,1,. ..,0) 



{1,1,..., 1). (4) 



This never happens if p 1. Things are thus very different 
when this fully symmetric scenario is assumed. 

It can be shown that in this fully symmetric scenario the 
competitive exclusion principle applies. This principle states 
that there can not coexist more populations than different re- 
sources (or ecological niches) in the long term if these pop- 
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R S 

Figure 10: (a) Mean absorption time T0 (in number of invasions) starting from the empty community 0, and mean number of species S as a function of the resource 
saturation R, showing a roughly linear growth for both of them, (b) Mean absorption time vs. 5 . 



ulations depend linearly on the resources ^Hofbauer and Sig- 



mund 1998 1. We can put this statement in mathematical terms. 



For the sake of simplicity, let us assume that there is a single 
trophic level with S species predating on the resource (at rates 
y+i, i - 1, ... ,5) and let us set a non-uniform direct compe- 
tition pij between pairs of species in that level. Let n, be the 
population density if species /, a,- its death rate in the absence 
of consumption and «o the amount of resource. The Lotka- 
Volterra equations for this system are 



(5) 



If the competition matrix is singular, we can find a non-trivial 
solution (ci,...,cs) for the linear system Yji^iPij = 0, y = 
I, . . .,S (note that, in particular, the fully symmetric scenario 
p = 1 renders the competition matrix singular). Hence 



^ c,(log n,) = ^ ciiy^itio - a,) = 

1=1 1=1 



(6) 



where we can assume that a is positive (otherwise change the 
sign of the c,). Integrating from to f we obtain 



Y[ni{tr-^Ce 



(7) 



This means that one of the densities must vanish in the limit 
f ^ 00, which proves competitive exclusion. 

There is a peculiarity of our model, though. If p = 1 the pop- 
ulation of the invader at equilibrium will not be zero because 
in our model all constants are uniform, so the equation to solve 
for c, is 2, c; = 0. This yields a = and spoils the argument. 
However, we have shown that, with our procedure of species 
extinction, the invader's population ends up below iic hence not 
being viable. This restores competitive exclusion, albeit in a 
weaker sense. The result Q is just a manifestation of this fact. 

It is important to notice that, for a non-singular competition 
matrix, the competitive exclusion principle is not guaranteed to 
hold. In particular, if p < 1 the intra- and interspecific compe- 
tition will have different magnitude, and the matrix of elements 



p„ - 1 and Pij - p {i j) will be non-singular. The argu- 
ment above does not apply anymore and, as a matter of fact, by 
integrating the equations for population dynamics we actually 
obtain more than one species coexisting with a single resource 
in the system. 

The interesting point brought about by the above discussion 
is that interspecific competition induces de facto a niche sep- 
aration for the species of the same level — which are therefore 
competing for the same resources — that allows them to circum- 
vent the competitive exclusion principle [for a more thorough 
discussion of this point see |Bastolla et al.|p005a|b[ )]. 

3.4. Absorption times 

So far we have discussed properties of the recurrent set of the 
Markov chain associated to the assembly process. But we have 
not discussed the possibility that the process may keep trapped 
for a long time in transient states. In order to ascertain this, we 
have calculated the mean absorption time from the empty com- 
munity to the end state. This can be done given the structure 
that the transition matrix P acquires after a permutation that 
reorders recurrent and transient states. We can thus write the 



matrix in a block form ( |Karlin and Taylor 1975 1 



P = 



U 
W V 



(8) 



where matrix U contains the transition probabilities within the 
recurrent set, and V contains transition probabilities between 
transient states. The average time from the transient state 
to reaching for the first time the recurrent set at state j is the 
element f,j of matrix T, where 



T = J^«y"-iw = ii-vy 



-w. 



(9) 



I being the identity matrix. This expression counts as n the 
absorption time when the process remains n - 1 time steps 
within the transient subset and jumps to a recurrent state in 
the n-th step. The mean absorption time for a process start- 
ing from the transient state / will thus be t, = YijeK Uj - (Ju),, 
where u - (1, . . . , 1)"'^. Since P is stochastic, 2j(V,7 + Wy) = 




Figure 11: (a) Probability n(m) that an invasion causes the extinction of at least a fraction m of the species of the invaded community, for several values of the 
resource saturation R. There is a characteristic magnitude j8"' around 1%. In order to have enough statistics, we have chosen values of R within the region where 
the number of communities in the recurrent set is above 1000 [see Figure 3 in Capitan et al. (2009 1]. Thus we have better statistics to compute the histograms than 
for smaller R. (b) The same probability but for transient states, nt(nj), for two values of R. The characteristic size of the avalanches here is about 2%. 



(yu)i + (Wu)i = 1 for all i e g - so Wu ^ (I - V)u or, 
equivalently, (I - Vy^Wu - u. Together with (|9| this implies 
(I - V)t - u, so solving this sparse linear system yields the 
absorption times for any transient state. Note that these times 
are proportional to ^ because of the form ([T} of our transition 
matrix. 



In Figure lOi we plot the mean absorption time t^i to reach 



the recurrent set starting from the empty community, along with 
the mean number of species S , which measures the size of the 
system. Both of them grow almost linearly with R, hence is 
roughly linear with S as well (see Figure [T0|^). Since the num- 
ber of states in each chain grows as e" the number of states 
of the Markov chain is very large compared to ^T0. Therefore 
the mean time to the end state is small compared to the system 
size. 

This result should be taken with a grain of salt, because it 
strongly relies on our assignment of probabilities to transitions. 
This, in turn, assumes that there is always availability of in- 
vaders, which may not be true if invaders come from a finite 
pool. The lack of potential invaders when the community is al- 
most "full" would decrease the probability of a new invasion 
and accordingly would increase the time that the process needs 



to reach the end state. What the result of Figure 10 1 is actually 



telling us is that the assembly graph is dominated by pathways 
in which most invasions are accepted. 

3.5. Extinctions distribution 

As we have previously described, the assembly process can 
be regarded as if the ecosystem self-organizes into a state re- 
sistant to invasions. Either for transient or recurrent states, the 
community is continuously undergoing avalanches of extinc- 



tions caused by new colonizations. Figure [TT^ shows a statis- 
tics of such avalanches in some recurrent sets. It represents the 
probability n(in) that an invasion causes an avalanche of mag- 
nitude greater than m (understood as the fraction of species that 
go extinct), averaged over the stationary state. We can see in the 
figure that this probability shows an exponential decay, with a 
typical avalanche size yS"^ of about 1% of the community, /3 be- 
ing the slope of the distributions in log-linear scale. Invasions 
never cause big perturbations in the community. 

We can calculate a similar distribution for the avalanches of 
extinctions in the transient states. Now we have to weight the 
magnitudes with the average fraction of visits to each transient 
state. Let us denote as Zij the average number of visits to state j 
starting from state /. The matrix Z - (zij) is then given by 



z = y" = (I - y)- 



(10) 



n=0 



Thus the number of visits to the transient /' starting from is 
- (u0Z)j, U0 being the row vector 6,0 (with a number of en- 
tries equal to the number of transient states). We can calculate 
( by solving the linear system 



^(i-y) = M0. 



(11) 



The resulting probability Tlt{m) that an invasion causes the 
extinction of at least a fraction m of the species of the invaded 
transient community is showed in Figure [TT]^. We also find 
an exponential behavior for the cumulative distribution, in this 
case with a mean characteristic fraction of species loss of 2% 
for transient avalanches. The species loss caused by invasions 
in the transient part of the graph is always small. 
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Figure 12: (a) Probability of invasion vs. mean number of invasions f) for R = 430 (witii a complex end state). Inset: probability of reconfiguration after invasion 
vs. mean number of invasions, (b) The same as (a) but for R = 540 (the end state is a single community), (c) Probability of invasion (dashed lines) and average 
number of species loss (full lines) vs. mean number of invasions, for two values of R with complex end states. At the time these two magnitudes coincide, a 
stationary state (the recurrent set) is reached. 



3.6. Time averages 

Computing the time evolution of averages is very simple, 
given the transition matrix P and some initial probability distri- 
bution (KarUn and Taylor 1975 1 — which in our case is simply 



the vector U0, since the assembly process starts from the empty 
community. We just need to calculate the power P' to obtain 
the transition probability matrix after t time steps. Thus we can 
obtain the probability of rejecting the invader at discrete time f 
as 

^r(f) = X^^v(n'^> (12) 



and that of accepting the invader as 



(13) 



\k-j\=l 



where the inner sum runs over transitions from j in which the 
invader is accepted. Obviously, the probability that the com- 
munity undergoes a reconfiguration because of the invasion is 
obtained as PJt) - 1 - P,(t) - Pi(t)- Figures [T2|i and \T2\ ) 
represent the dependence in time of the probabilities Pi and P.^ 
in two cases: one with a a complex end state (a), and another 
with a single community as end state (b). Notice that all curves 
collapse, for small ^, when divided by ^ and plotted against 
(mean number of invasions). 



In Figure 12 ; we show the probability of invasion f and 



the average species loss defined as 



(14) 



where (A5),y is the species loss in the transition from j to k 
and the prime denotes that we ignore in the sum transitions in 



which the invader is accepted. When these two magnitudes are 
equal there is an equilibrium between the average frequency of 
invasions and the average number of species loss. This is a fin- 
gerprint of the reaching of the stationary state. As expected, this 



time is comparable to the absorption time shown in Figure 10 1. 



Another important magnitude is biodiversity. Figure 13 i rep- 



resents the evolution of the average number of species for sev- 
eral values of R. In all cases, this average number grows mono- 
tonically until reaching the stationary state, so biodiversity and 
resistance to invasion are positively correlated, in agreement 



with previous assembly models ( Law and Morton 1996 Mor- 
|tonandLaw||T997l ). 

Figure \T3[ > represents the evolution of the total population 
density B of each community. If we assume, for the sake of 
simplicity, the same weight per individual for all species in our 
model communities, then B can be regarded as the total biomass 
in the community. Although there is a clear trend for biomass to 
increase, it is not always at its optimum in the stationary state. 
This is very clear in the figure for R = 470, a value at the onset 
of the appearance of the fourth trophic level. This agrees with 



the analysis performed by Virgo et al. ( 2006 1 on their assembly 
model. 

We have also studied the time dependence of the average 
number of trophic levels during the assembly, which is shown 
in Figure [13};. At R = 470 the process stays a certain time 
trapped in three-level communities until the fourth level is fi- 
nally accepted. This effect becomes lower upon increasing R, 
until there is no trapping and the fourth level is reached straight 
away. 



Figure 131 shows a typical time evolution of the average re- 



turn time along the assembly until reaching the stationary state. 
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Figure 13: Species richness (a) and total population (b) vs. mean number of invasions for several values of the resource saturation S. At S = 470 (showed with 
dash dotted lines) the ecosystem crosses over from 3 levels to 4 levels (this crossover corresponds to the non-monotonic behavior of the total population), (c) Mean 
number of levels vs. mean number of invasions. At the onset of the fourth level, the ecosystem stays some time trapped in three-level communities, (d) Typical 
variation of the average return time with f /. 



Communities are less resilient (have larger return time to equi- 
librium) as time increases. Thus, there is a trade-off between 
robustness (resistance against invasions) of the ecosystem and 
dynamic stability which is resolved by sacrificing the latter in 
favor of the former. 

4. Discussion 

In this work we have provided a full account of results of 



the model introduced in |Capitan et al. (2009j). The results pre- 
sented here have been obtained from a direct analysis of the 
Markov chain describing the assembly process. Among the 
novel results here presented are the dependence of many bio- 
logical observables on the amount of available resource (imple- 
mented through parameter R), the average times that the process 
needs to reach the recurrent set, the statistics of avalanches in 
both transient and recurrent states, and the time evolution of any 
observable. 

Our model might be considered as a benchmark of the assem- 
bly process that builds up ecological communities. As such, we 
do not aim at providing a realistic description of an ecosystem 
but at capturing, in a very simplified model, the essential mech- 
anisms that do occur in the construction of real ecosystems. The 
model rests on some oversimplistic features: communities are 
strictly organized in levels, predation occurs only between con- 
tiguous levels, species of a given level are trophically equiva- 
lent, model parameters are chosen uniformly and the popula- 



tion dynamics is ruled by simple Lotka-Volterra equations. In 
spite of this, our model exhibits the same behavior as all other 
assembly models reported in the literature. This indicates that 
this behavior is very robust, and probably shared by real sys- 
tems and simple models alike. 

Thanks to these oversimplifications the model provides im- 
portant advantages on previous assembly models. The main 
one is that we can trace all pathways of the assembly process. 
This allows us to compute exactly all the observables of a com- 
munity and to characterize in a very precise manner the sta- 
tionary state of the ecosystem. Our model also has a species 
pool, as standard assembly models, but because we allow for 
an arbitrary number of trophically equivalent species, the pool 
is infinite and the model does not suffer from the problem of 
exhaustion of good invaders that may trap the community in a 
transient state (I Case| 1 1 99 1 [ [Levine and D' Antonio] 1 1 999) 1 . This 
has permitted us to build communities with hundreds of species 
and explore the influence of different elements on the behavior 
of the assembly process. 

Therefore, we are not limited, as in standard assembly mod- 
els, to compute averages over a set of realizations of the pro- 



cess. As we pointed out in Capitan et al. (2009) , the number 
of shortest pathways leading from to the recurrent set can be 
enormous. For instance, for R - 300 (a case with an absorb- 
ing community of three trophic levels and 50 species), there are 
~ lO"' different minimum-length pathways. This number is far 
from anything a simulation can come close to. 
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There is, of course, a concern about having trophically equiv- 
alent — hence indistinguishable — species. The grouping of 
trophically equivalent species is a common practice in study- 
ing food webs, so it is tempting to do so in this model. If 
we do it, the model becomes equivalent to a chain, for which 



Lotka-Volterra dynamics is well characterized (Hofbauer and 



Sigmund 1998) , and the invasion process seems to become 
trivial. This is not true, though: if p + 1, i.e. if intra- and 
interspecific competition are different in magnitude, intraspe- 
cific competition in the equivalent chain explicitly depends on 
St, so invasions modify the parameters of the chain and the inva- 
sion process becomes non trivial. Thus, it is because of the di- 
rect interspecific competition p < 1 that this equivalence breaks 
down and the model departs from triviality. We have explicitly 
shown that choosing p = 1 brings about the competitive exclu- 
sion principle, and indeed the model turns into a chain. But for 
any p < 1 this does not longer hold. Interspecific competition 
is thus an effective way of creating new niches. 

Let us now summarize the main conclusions we can extract 
from the present analysis of the model. 

As our model ecosystems evolve we observe three trends: 
biodiversity increases, resistance to invasion increases and all 
species decrease their populations. In the steady state biodiver- 
sity is at its maximum, all populations are close to the extinction 
level and either invasions are rejected or they produce transi- 
tions between a set of communities with a very similar struc- 
ture. All three features are related. The increase in biodiver- 
sity is unavoidable because of the constant flux of colonizers; 
however, as the number of species increases, their populations 
necessarily decrease because all share the same resource. The 
invasion process guarantees that this is done in the most effi- 
cient way, because inefficient invasions cause extinctions in the 
community and force a more equilibrated rearrangement of the 
populations. This, in turn, justifies the increasing resistance to 
new invasions. At the end, all populations are so close to ex- 
tinction that either no new invasions are possible, or they just 
cause small rearrangements that leave the community in a sim- 
ilar state. 

Final communities have typically three or four trophic levels 
— only ecosystems with more than 200 species generate five 
trophic levels. On the other hand, the number of species in each 
level has a pyramidal structure. Both features are in qualita- 
tive agreement with what is observed in real ecosystems ( ,Cohen| 
et al. 1990| l and we have discussed at length the properties of 
the population dynamics equations that explain these features 
inp 



Capitan and Cuesta. ( .201O). 



As already advanced in Capitan et al. ( 2009| l the end state 
is always unique, and this is consistent with previous assembly 
models (Mort on and Law|[T997j . However, there is a caveat 
that should be made on this point related to the indistinguisha- 
bility of species within the same trophic level: the end state 
is unique as long as we consider only the number of species 
at each level. Whether two communities with the same num- 
bers have the same or different species is meaningless for this 
model, so the conclusion is not definitive. In fact, some rel- 
atively recent experiments on aquatic microbial communities 
establish that productivity-biodiversity relationships depend on 



the history of assembly (Fukami and Morin 2003| l, and it is 
our guess that the independence on history resulting from this 
model might be an artifact of the indistinguishability of species. 
Refined versions of this model may clarify this issue. 

As for the robustness of the above results, we have tried other 
values of the direct competition parameter, namely p = and 
p = 0.7, to test its influence. No qualitative difference with the 
behavior reported here is found. Nonetheless, there are three 
quantitative effects that we have observed as p increases: resis- 
tance to invasion increases, appearance of new trophic levels is 
hindered and the number of communities in complex end states 
decreases. Varying j_ has similar effects; in fact, the product 
y+y- - O.ly^ provides a quantitative estimate of indirect com- 
petition. 

It can be argued that parameters should depend on the trophic 
level rather than being uniform for all species. It is very easy 
to show that this does not change the dynamic stability patterns 
because in that case one can also construct a Lyapunov function 



[see [Capitan and Cuesta (2010 1]. We have not attempted any 
test in this respect, but it is hard to believe that such a variant 
of the model will produce any qualitative difference. The as- 
sembly graphs will be similar to the ones found for the present 
model. Something more can be said about the invasion rate. We 
have presently assumed that the invasion probability is the same 
for all trophic levels, but notice that the assembly graph is ut- 
terly independent on this choice, so certainly choosing a differ- 
ent invasion probability will change the numerical value of the 
nonzero entries of the transition matrix P, but only them. The 
graph, as well as the structure of transient and recurrent states 
of a finite Markov chain, only depends on which elements of 



P are zero (KarUn and Taylor 1975 1, so not just the graph but 
the set of communities in the end state will be exactly the same 
as those reported here (the probability distribution in the steady 
state will, of course, be different). 

Perhaps the most important limitation of this model is the 
choice of the Lotka-Volterra equations. The choice of popula- 
tion dynamics has been reported to have a strong influence in 



the final shape of ecological communities (Drossel et al. 2004 



Lewis and Law[[2007[ ). Introducing non-linear equations leads 



to more complex stability patterns than simply rest points. How 
to account for them is not yet clear to us, but neither is whether 
this will really affect the qualitative behavior of the assembly 
process. Thus, this remains an important open question that 
deserves further analysis. 

There are further open questions such as the application of 
this model to metacommunities. The resulting ecosystems can 
be readily altered when migration takes place among spatially 
distributed patches. With a simple model like ours, it might 
be possible to build up an assembly graph between different 
communities in different patches. The interplay between com- 
munities in different patches could lead to an outcome different 
from the one we obtain with a single patch. On the other hand, a 
simple model like this can provide us with basic understanding 
of complex processes such as, for instance, the rebuilding of a 
natural community after its degradation. Very little is known 
about the processes that helps to reconstruct damaged commu- 
nities, and a simple framework like ours could provide some 
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hints about how to tackle this problem from a theoretical point 
of view. 

The final take-home message from this work is this: we 
should not be afraid of oversimplifications in complex systems. 
Complexity normally arises as a consequence of a collective 
behavior of many entities, not as a result of the complexity 
of interactions. The key point is whether we are retaining the 
basic ingredients yielding the desired output. We have shown 
that there is no qualitative difference between the results of this 
oversimplified model and previous, more sophisticated assem- 
bly models. And there is a lot to gain from the wider view that 
this model provides of the process and the much higher control 
we have on the parameters. Many questions that are hard (or 
even impossible) to answer in previous model have a clear-cut 
answer here. And even if they may be too simplistic, they can 
still guide our intuition when dealing with real ecosystems. 
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